#ifndef ADAPTIVE_PRECISION_H
#define ADAPTIVE_PRESICION_H

#include <utility>

std::pair <double, double> sum(double a, double b)
{
   double x = a + b;
   double b_virtual = x - a;
   double a_virtual = x - b_virtual;
   double b_roundoff = b -b_virtual;
   double a_roundoff = a - a_virtual;
   double y = a_roundoff + b_roundoff;
   return std::pair<double, double>(x, y);
}

template <size_t s>
std::pair<double,double> split(double a)
{
   double c = ((1LL << s) + 1LL) * a;
   double ab = c - a;
   return std::pair<double, double>(c - ab, a - c + ab);
}

template <size_t m>
void grow_expansion(double* e, double b, double* h)
{
    for(size_t i = 0; i < m; i++)
    {
       std::pair<double,double> temp = sum(e[i], b);
       b = temp.first;
       h[i] = temp.second;
    }
    h[m] = b;
}

template <size_t m, size_t n>
void expansion_sum(double* e, double* f)
{
    for(size_t i = 0; i < n; i++)
    {
        grow_expansion<m>(e + i, f[i], e + i);
    }
}

std::pair<double,double> mul(double a, double b)
{
    double x = a * b;
    double a_hi, a_lo, b_hi, b_lo;
    std::pair<double, double> a_int = split<(std::numeric_limits<double>::digits - std::numeric_limits<double>::digits / 2)>(a);
    a_hi = a_int.first;
    a_lo = a_int.second;
    std::pair<double, double> b_int = split<(std::numeric_limits<double>::digits - std::numeric_limits<double>::digits / 2)>(b);
    b_hi = b_int.first;
    b_lo = b_int.second;
    double e1 = x - (a_hi * b_hi);
    double e2 = e1 - (a_lo * b_hi);
    double e3 = e2 - (b_lo * a_hi);
    double y = (a_lo * b_lo) - e3;
    return std::pair<double, double>(x, y);
}

template <size_t n>
int get_sign(double * e)
{
    for (int i = n - 1; i >= 0; i--)
    {
        if (e[i] > 0)
        {
            return 1;
        }

        if (e[i] < 0)
        {
            return -1;
        }
    }
    return 0;
}
#endif